﻿      subroutine polarlongtide(t,pole)
!Tidal variations in polar motion and polar motion excitation due to
!long period ocean tides
      implicit none
	integer::i,j
	real*8::t,pi,RAD,th,pole(4),f(5)
	real*8::frk(10,14),hp,hr,fp,fr
      data(frk(i,1),frk(i,2),frk(i,3),frk(i,4),frk(i,5),frk(i,6),frk(i,7),frk(i,8),
     * frk(i,9),frk(i,10),frk(i,11),frk(i,12),frk(i,13),frk(i,14),i=1,10)/
     *  1,0,2, 0,1,    9.12,  4.43,-112.62,  5.57, 21.33, 205.83, 67.21, 269.95, 21.17,
     *  1,0,2, 0,2,    9.13, 10.72,-112.56, 13.48, 21.30, 497.59, 67.27, 652.59, 21.14,
     *  0,0,2, 0,1,   13.63, 27.35, -91.42, 30.59, 13.31, 841.32, 88.42,1002.12, 13.15,
     *  0,0,2, 0,2,   13.66, 66.09, -91.31, 73.86, 13.27,2028.73, 88.53,2414.94, 13.11,
     *  0,0,0, 2,0,   14.77,  5.94, -87.13,  6.42, 11.75, 168.13, 92.70, 194.74, 11.60,
     *  1,0,0, 0,0,   27.56, 43.74, -56.70, 31.12, -0.91, 643.61,123.13, 520.16, -1.06,
     * -1,0,0, 2,0,   31.81,  8.85, -51.11,  5.42, -4.21, 111.62,128.72,  79.23, -4.36,
     *  0,0,2,-2,2,  182.62, 86.48, -20.30, 99.77,175.57, 118.56,159.42, 336.32,175.46,
     *  0,1,0, 0,0,  365.26, 17.96, -17.38,152.15,170.60,   3.33,161.60, 332.53,170.51,
     *  0,0,0, 0,1,-6798.38,208.17, 166.89,186.98,166.67, 221.43,166.88, 175.07,166.68/
!---------------------------------------------------------------------------
   	pi=datan(1.d0)*4.d0;RAD=pi/180.d0;pole=0.d0
      !Delaunay variables l, lp, F, D, omga (IERS2010, 5.43)
	f(1)=134.96340251d0+(1717915923.2178d0*t+31.8792d0*t**2
     >	     +0.051635d0*t**3-0.0002447d0*t**4)/36.d2
	f(2)=357.52910918d0+(129596581.0481d0*t-0.5532d0*t**2
     >	     +0.000136d0*t**3-0.00001149d0*t**4)/36.d2
	f(3)=93.27209062d0+(1739527262.8478d0*t-12.7512d0*t**2
     >	     +0.001037d0*t**3-0.00000417d0*t**4)/36.d2
	f(4)=297.85019547d0+(1602961601.209d0*t-6.3706d0*t**2
     >	     +0.006593d0*t**3-0.00003196d0*t**4)/36.d2
	f(5)=125.04455501d0+(-6962890.5431d0*t+7.4722d0*t**2
     >     +0.007702d0*t**3-0.00005939d0*t**4)/36.d2
	f=dmod(f*RAD,2.d0*pi)
	do i=1,10!m=0长周期10
	  th=0.d0
	  do j=1,5
	    th=th+f(j)*frk(i,j)
        enddo
        th=-th
        hp=frk(i,7);fp=th+frk(i,8)*RAD;hr=frk(i,9);fr=-th+frk(i,10)*RAD
        pole(1)=pole(1)+hp*dcos(fp)+hr*dcos(fr)
        pole(2)=pole(2)-hp*dsin(fp)-hr*dsin(fr)
        hp=frk(i,11);fp=th+frk(i,12)*RAD;hr=frk(i,13);fr=-th+frk(i,14)*RAD
        pole(3)=pole(3)+hp*dcos(fp)+hr*dcos(fr)
        pole(4)=pole(4)+hp*dsin(fp)+hr*dsin(fr)
	enddo
	return
	end
